预测海藻数量案例
本章主要介绍了一些数据挖局的基本任务:数据预处理、搜索性数据分析、预测模型的构建
问题描述以及目标
河流当中有害藻类的含量会极大影响河流的生态环境,本案例的目的是为了通过使用过去一年采集的不同时间段的水样成分、季节、流速、河流大小这些数据建立模型来预测海藻数量
数据说明
案例数据来自ERUDIT,案例总共有两个数据集。
第一个数据集有200个水样,记录了同一条河流在概念同一季度的三个月内收集的水样的平均值,每天记录由11个变量构成,其中3个名义变量它们描述了水样采集的季节、河流大小、流速,余下的8个变量是不同的化学物质含量。
- 最大pH值
- 最小含氧量(O~2~)
- 平均氯化物含量(CL)
- 平均硝酸含量(NO~3~^-^)
- 平均氨含量(NH~4~^-^)
- 平均正磷酸盐含量(PO~4~^3-^)
- 平均磷酸盐含量(PO~4~)
- 平均叶绿素含量
第二个数据集来自140个额外观测值。他们的基本结构和第一个数据集一样,但是它不包含7种藻类的频率数目,这些额外的观测值可以视为测试集,案例的主要目标是预测140个水样中7种藻类的频率,所以这是一个预测性数据挖掘类型任务,在这类问题中任务时间里预测模型,并预测在给定预测变量取值是相应的目标变量的值。预测模型也可能会说明哪一个预测变量对目标变量有较大的影响,即模型可能提供影响目标变量因素的一个综合描述。
数据加载到R
我们考虑两种数据载入R的方法,第一种是利用教材提供的R添加包直接加载另外一种方案是下载数据的文本文件之后将文件读入内存。
直接使用lib加载
1 | > library(DMwR) |
我们通过head看到了数据的前6行
读取文本文件
以文件名A.txt为例,文件包含200行、每列使用空格切分,缺失数值用字符串‘XXXXXXX’来表示
1 | algae <- read.table('A.text', header = F,dec = '.',col.names = c('season','size','speed','mxPH','mnO2','Cl','NO3','NH4','oPO4','PO4','Chla','a1','a2','a3','a4','a5','a6','a7'),na.strings = c('XXXXXXX') |
header参数表示要读取的文件第一行不包含变量名(列名),dec表示在数值位当中用来区分小数位的字符,col.name为正在读取的变量提供一个名称向量来表示变量名(因为文件第一行没有提供),最后na.string设置默认字符串
数据可视化和摘要
R当中提供了summary函数来获取统计特性的一个摘要信息
1 | > summary(algae) |
从结果中我们可以看到在概览展示中,对于名义变量(因子)他给出了每个可能取值的频次(season字段),而对于数值型变量则直接提供了均值、中位数、四分位、极值、和空值的数量 例如Cl数据总共有10个缺失。
下面我们绘制一下mxPH的直方图
1 | hist(algae$mxPH,prob = T) |
设置参数prob=T表示我们要显示的是概率,否则会显示为数量,通过图形我们可以看到mxPH的分布非常接近正态分布,它的值都聚集在均值周围,通过Q-Q图来检测该变量是否为正态分布,下面需要用到的是R的添加包car当中所提供的函数qq.plot
1 | library(car) |
最终结果如下图所示
使用箱图进行数据检查
使用统计图进行数据监测可以使用箱图来进行,如下代码:
1 | #绘制箱型图 |
最终会结果如下图
箱图包含的信息很多,能够非常直观的看出来数据的范围、分位数、平均数、中位数、上边缘、下边缘、异常值,wiki链接
有时候当有离群值时,可以通过下面这种方式来进行观测
1 | #绘制数据点 |
我们会得到如下图线
我们可以直接使用identify来将点击采集的数据点保存下来
1 | plot(algae$NH4,xlab = '') |
还可以使用添加包lattice来进行数据检查,该包提供了大量优秀的图形工具。
条件绘图是依赖于某个特定因子的图形表示,因子是一个取值为有限集合的名义变量,例如,对于变量size的不同区直,可以绘制变量a1的一组箱图,每个箱图对应于变量size的某个特定值的水样子集,通过这些箱图可以研究size变量是如何影响变量a1的分布的。
1 | library(lattice) |
绘制图形如下
类似的还有分位箱图,绘制代码如下
1 | > library(Hmisc) |
分位箱图当中的点代表了不同大小的河流中海藻频数的均值,图中的竖线分别代表了变量的第一分位数,中位数,第三个分位数,小竖线则代表了数据的实际取值,这些值的分部信息则由分位数图来体现,该图提供的信息是要多余传统箱型图的,我们可以确认小型河流有更高频率的海藻,但是我们也观察到小型河流的海藻频率的分布比其他类型河流的海藻分布要分散
这种类型的条件绘图不局限与名义变量也不局限于单个因子,只要先把连续变量“离散化”也同样可以进行条件绘图,下面给出一个两个因子的条件绘图的例子。考虑a3在给定变量season和变量mnO2下的条件绘图,变量mnO2是一个连续变量,绘制代码如下:
1 | # 将变量mnO2离散化为因子类型,number设置区间个数,overlap设置两个区间之间靠近边界的重合度,na.omit提出空值 |
数据缺失
通过观察数据我们可以看到有一些变量是有缺失值的,这种情形在现实问题中非常普遍,这会导致一些不能处理缺失值的分析方法无法应用,以下是一些处理缺失值数据的几种常用策略:
- 将含有缺失值的案例剔除
- 根据变量之间的相互关系填补缺失值
- 根据案例之间的相似性填补缺失值
- 使用能够处理缺失值数据的工具
最后一种方式最为严格,它直接限制了我们能够使用的工具。
将缺失部分剔除
剔除含有缺失数据的记录非常容易实现,尤其是当这些记录所占的比例在可用数据集中非常小的时候,这个选择就会比较合理。
首先我们获取包含有缺失值的数据个数:
1 | # complete.cases产生一个布尔值向量,该向量的元素个数与algae数据中的行数相同,如果该行数据不含有NA值(完整数据)函数返回值为T |
通过输出我们可以看到62和199行记录中11个解释变量有6个是缺失值,在这种情况下我们可以直接将其进行剔除
1 | # 使用c生成一个向量取负数表示丢弃 |
我们还可以计算数据当中每行的空数据个数
1 | apply(algae,1,function(x) sum(is.na(x))) |
apply被称为元函数,他们会把后面传递的函数传递给数据,据此我们可以自定义函数将数据缺失达到一定比例的数据剔除
1 | function (data, nORp = 0.2) |
利用高频率值来填补缺失值
填补含有缺失值记录的另一个方法是尝试找到这些缺失值最可能的取值。填补缺失数据最简单和快捷的方法是使用一些代表中心趋势的值,代表中心趋势的值反映了变量分布的最常见取值,因此中心趋势值是最自然的选择。有多个代表数据中心趋势的指标例如:平均值、中位数、众数。最合适的选择由变量的分布决定,对于接近正态分布的值来说,所有数值都较好的分布在平均值周围,平均值就是最佳选择,然而对于偏态分布或者含有离群值的变量来说,选择平均值就不好了。偏态分布的大部分值聚集在平均值分布的一次,因此平均值不能作为最常见值的代表,另一方面,离群值的存在也影响了平均值。这就导致了平均值不具有代表性,因此我们在进行操作之前应该首先观察数据分布。对于偏态分布或者有离群值的分布而言,中位数就是一个更好一些的代表数据中心趋势的指标。
例如在样本algae[48,]当中mxPH就有缺失值,由于该变量分布近似正态分布,我们可以使用平均值来填补这个洞。
1 | algae[48,'mxPH'] <- mean(algae$mxPH,na.rm = T) |
而对于变量Chla这个变量在12行上有缺失值。这也是平均值不能代表大多数变量值的一种情况,事实上Chla的分布偏向于较低的数值,并且他有几个极端值,这些使得平均值(13.971)不能代表大多数的变量值。因此我们使用中位数来填补这一类的缺失值:
1 | algae[is.na(algae$Chila),'Chla'] <- median(algae$Chla,na.rm = T) |
由于缺失值的存在会导致某些方法不能使用,所以使用上面的方法填补缺失值通常也认为不是很好的方法。虽然以上方法速度快,适用于大数据集,但是他可能导致较大的数据偏差,影响后期的数据分析工作。然而使用五篇方法寻找最佳数据填补值复杂,对于大型数据挖掘问题可能并不适用。
通过变量的相关关系来填补缺失值
另外一种获取缺失值较少偏差估计值的方法是探寻变量之间的相互关系。比如,通过变量值之间的相关关系,能够发现某变量与mxPH高度相关。这就使得我们得到更好更可能的填补值来填充到底48条记录当中。
使用如下命令获取变量间的相关性
1 | cor(algae[,4:18],use = 'complete.obs') |
函数cor功能是产生变量之间的相关值矩阵,之所以从第四个变量开始是因为前三个是名义变量(season/size/speed),设定参数complete.obs在计算时会忽略掉含有NA的记录,相关值在1或者-1周围表示相应的两个变量之间有强正或者强负的线性相关关系。之后可以在其它函数当中得到变量间线性相关的近似函数形势。
函数cor的输出并不是很清晰,我们可以通过函数symnum来改善结果的输出形式
1 | symnum(cor(algae[,4:18],use = 'complete.obs')) |
通过谷歌文档找到了对于symnum输出的解释,图中的matrix被填上了各个字符,注意他不是空哪个是空白字符,整张图表示0~0.3之间的数值使用空格表示,0.3~0.6之间数值使用’.’表示,以此类推,从图中我们可以看到oPO4和PO4是强相关的大于0.9,因此对于PO4和oPO4我们可以利用他们的相关性来填补他们的缺失值,首先我们要找到这两个变量的线性相关关系。
1 | > algae_filter <- algae[-manyNAs(algae),] |
使用函数lm()可以获得形如Y = β~0~ + β~1~X~1~ + ··· + β~n~X~n~ 的线性模型。在本例中就是PO4 = 42.897 + 1.293 x oPO4,利用此变量可以对数据进行补充,我们可以通过自定义一个函数来使用sapply来对数据进行填充。
1 | algae_new <- algae[-manyNAs(algae),] |
sapply第一个参数是一个向量,第二个参数是一个函数,sapply会使用函数参数遍历向量并最终生成一个新的向量返回。
对于线性关系的研究使得我们能够填充一些新的缺失值,然而还有几个观测值含有缺失值,可以试着探索案例数据中含有缺失值的变量和名义变量之间的关系,可以通过R添加包lattice中的函数来绘制条件直方图来进行
1 | histogram(~mxPH | season,data = algae) |
代码绘制了不同季节变量mxPH的直方图,每个直方图对应于某个季节的观测值数据,季节不是按照自然地时间顺序,因为没有给定绘制变量的顺序,我们可以通过如下来进行指定
1 | algae$season <- factor(algae$season,levels = c('spring','summer','autumn','winter')) |
我们可以通过直方图拓展到多个变量来进行观察
1 | histogram(~mxPH| size *speed,data =algae) |
图形绘制了河流大小和速度的组合对mxPH值的变化,需要说明的是它没有说明速度较低切规模较小的河流的相关信息,因为具备这个信息的恰恰是丢失了信息的第48个样本
通过另一种方式也可以得到类似的信息,但这次应用变量的具体取值如下:
1 | stripplot(size ~mxPH|speed,data = algae,jitter = T) |
结果如下图所示,jitter = T表明对Y轴的变量进行小范围的随机排列,这样可以避免由于相同值之间相互重叠而导致失去某些特定值的集中度信息。
这种类型的分析可以用到其他含有缺失值的变量中,而且这种分析会需要用到多种组合来进行分析,这种方法可以用到少量名义变量的较小数据的分析当中
通过探索案例之间的相似性来填补缺失值
待续。。。
转载请注明来源链接 http://just4fun.im/2018/04/14/learning R(二)//) 尊重知识,谢谢:)